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Abstract 

Macroscopically populated quantum superpositions pose a question to what 
extent macroscopic world obeys quantum mechanical laws. Recently such super- 
' f-H '■ positions for light, generated by optimal quantum doner, were demonstrated. 

O | They arc of fundamental and technological interest. We present numerical meth- 

ods useful for modeling of these states. Their properties are governed by Gaus- 
sian hypergeomctric function, which cannot be reduced to neither elementary 
nor easily tractable functions. We discuss the method of efficient computation 
of this function for half integer parameters and moderate value of its argument. 
We show how to dynamically estimate a cutoff for infinite sums involving this 
function performed over its parameters. Our algorithm exceeds double preci- 
sion and is parallelizable. Depending on the experimental parameters it chooses 
one of the several ways of summation to achieve the best efficiency. Methods 
presented here can be adjusted for analysis of similar experimental schemes. 

Keywords: macroscopic quantum superpositions, macroscopic entanglement, 
optimal quantum cloning, Gaussian hypergeometric function, quantum optics 
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Programming language: C with OpenMP extensions (main numerical program), 
Python (helper scripts) 

Computer: modern PC (tested on AMD and Intel processors), HP BL2x220 
Operating system: Unix/Linux 

Has the code been vectorized or parallelized?: yes (OpenMP) 
RAM: 200 MB for single run for 1000 x 1000 tile 

Running time: l-2h for 1000 x 1000 tile, depending on the values of parameters 
Classification: 4.15, 18 

Nature of problem: Recently macroscopically populated quantum superpositions 
for light, generated by optimal quantum doner, were demonstrated. They are 
of fundamental and technological interest. Their properties are governed by 
Gaussian hypergeometric function 2-F1 of half-integer parameters, which cannot 
be reduced to neither elementary nor easily tractable functions. Computation of 
photon number distribution, visibility and mean number of photons, necessary 
for characterization of these states, requires evaluation of infinite sums involving 
this function performed over its parameters. 

Solution method: MQSVIS program suite computes various quantum indicators, 
such as photon number distribution, visibility, mean number of photons, vari- 
ance for macroscopic quantum superpositions of light. It takes losses (modeled 
with a beamsplitter) and imperfect photodetection (modeled with a Weierstrass 
transform applied to the photon number distribution) into account. Cutoffs of 
the infinite hypegeometric sums are estimated dynamically, and precision is en- 
hanced with computation of expressions in logarithmic form. Depending on 
the experimental parameters, program chooses one of the several ways of sum- 
mation to achieve the best efficiency. Program is parallelized using OpenMP 
standard, which ensures best utilization of multicore processors, and splits the 
work into tiles computed with different nodes of a computer cluster. This allows 
to compute the required indicators for realistic values of the parameters. 

2. Introduction 

Macroscopic quantum superposition (MQS) and entanglement were pointed 
out by E. Schrodinger in 1935. He posed a gedanken experiment where the 
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quantum formalism was applied to a macroscopic object - a cat and a micro- 
scopic object - a two-level radioactive atom, describing them in a joint quantum 
superposition [l[ . These objects were closed in a box and if the atom decayed 
at a random time, additional mechanism installed inside was activated to kill 
the cat. The conclusion was striking: the cat was dead and alive at the same 
time, as long as no one opened the box. It revealed the phenomenon of quantum 
entanglement: the state of the cat is random (dead or alive) but perfectly cor- 
related with the state of the atom (decayed or not). It also planted a question 
to what extent macroscopic objects obey laws of quantum mechanics. 

Small MQS were produced in: nanoscale magnets, laser-cooled trapped ions, 
photons in a microwave cavity, Ceo molecules, superconducting devices [2 and 
macroscopic-size diamonds cooled 0] and in room temperature (2011) [4|. In 
2007, macroscopically populated polarization entangled superpositions of light 
were generated in the room conditions by optimal quantum cloning (OQC). It 
is based on parametric frequency down conversion (PDC) di where a higher 
energetic blue photon in a given polarization is turned, with some probability, 
into two lower energetic red photons by a nonlinear crystal preserving the polar- 
ization. The more the crystal is pumped with the blue photon beam, the higher 
the probability of conversion, resulting in multi-photon output. OQC does not 
violate the no-cloning theorem Q , which states that an unknown quantum state 
(here: polarization) cannot be copied perfectly. It makes imperfect copies, char- 
acterized by a cloning fidelity lower than one There are two kinds of MQS 
of light produced by OQC: the micro-macro singlet state, where a single pho- 
ton is entangled with a "macroscopic qubit" , and the bright entangled squeezed 
vacuum 0, They contain 10 5 -10 13 photons. 

MQS are promising alternative for quantum technologies. They allow explor- 



ing the quantum-to-classical transition [lOl . 1 1 11 ] , efficient interaction with matter 
and photons (l2l |. studying the principle of quantum measurement and form 
a non-Gaussian class of quantum states, necessary for fault tolerant quantum 
computing fljj . Importantly, they give hope for a loophole- free Bell inequality 
test 



14| enabling the ultimate test of the quantum theory. 



Theoretical and experimental analysis of MQS is challenging due to their 
large complexity and fast decoherence both scaling exponentially with 

Hilbcrt space dimension. Commercial numerical tools (Matlab, Wolfram Math- 
ematica) and libraries (BLAS, NetLib, LAPACK) model MQS only for small 
populations. Intuitions based on these results are often misleading for high 
populations @, E^. Here, slowly convergent multiple infinite hypergeometric 
series arises, which is intractable for these applications. This results in lack of 
commonly accepted model of MQS. 

In this paper, we develop numerical model for MQS of light and discuss the 
method of efficient computation of Gaussian hypergeometric function for half- 
integer values of its parameters and moderate value of the argument. We show 
how to dynamically estimate a cutoff for infinite sums involving this function 
performed over its parameters. Our algorithm offers precision exceeding double 
precision by several orders of magnitude. Depending on the experimental pa- 
rameters it chooses one of the several ways of summation to achieve the best 
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Figure 1: (Color online) a) Phase-covariant optimal quantum cloncr. b) Poincarc polarization 
sphere. The blue great circle contains equatorial states of polarization. 



efficiency. It is parallclizable, which allows better utilization of multi-processor 
computers. We tested our numerical results with recently developed analytical 
model for MQS for certain regimes of parameter values (l4j . The numerical 
methods presented here can be adjusted for analysis of similar experimental 
schemes. 

This paper is organized as follows. In Section [3] wc introduce MQS, their 
properties and show the origins of hypergeometric function. In Section @] we 
discuss methods used in numerical model of MQS. Section [5] includes examples: 
computation of MQS photon number distribution and distinguishability. 

3. Theoretical background 

3.1. Entangled Macroscopic Quantum Superpositions 

In our analysis, we focus on the micro-macro singlet state, but our results 
can be generalized for the bright squeezed vacuum. The singlet is produced 
by phase covariant optimal quantum doner Q. This method requires first a 
pair of linearly polarized photons created in a usual state through parametric 
frequency down conversion (PDC) l/v / 2(|lff)a|ly)b — |lv)a|lff)&)- Subscripts 
a and b denote two spatial modes and H and V denote horizontal and vertical 
polarization, see Fig. OJ,. The equatorial states of the Poincare sphere of all 
polarization states, see Fig.Q}), are given by the following transformation of the 
linear polarization a\ = \/^2{e lv a ] H + e^a^), oL = i/V2{e lv a^ H - er^aJy), 

where aj, and a) ± are creators for two orthogonal polarizations tp and ip 1 - . This 
subspacc, parametrized by the polar angle ip € (0, 271"), is privileged for the phase 
covariant doner. Here, its Hamiltonian reads H = ihT ^(aj,) 2 + (a^) 2 ^) + h.c. 

and shows that all equatorial states are cloned equally well (the form of H is the 
same for all tp). Due to rotational invariance of the singlet, it can be expressed 
in this basis 1/V2(\l ip ) a \l tp ±}b — |l v ,-L) |l ¥J )h). Next, one of its spatial modes, 
e.g. b, is amplified by the doner to create a multi-photon state |$) = e lHt / h \l v ) 
or = e' lHt / h \l v ±) . This unitary evolution leads to the micro-macro singlet 

I*") = 1/V2(\l v ) a \$ ± ) b - | V) |$) 6 ). (1) 
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Figure 2: Probabilities 7? (solid line, left Y axis) and 7^ (dashed line, right Y axis) computed 
for gain g = 4. 



The multi-photon states are infinite superposition of photon number states 

00 oc 

l*>= 5>tf|(2i+l) v ,(2jV>, |$ x )= £ 7tf 1(2^,(2* + (2) 

i,j=0 i,j=0 

with the real- valued probability amplitude 

7 « = C- 2 (T 9 /2) i+ V(2i + l)!(2j)!/(i!j!) = C fl 2 7*070.,, (3) 

where ^^=0 7ij = 1- Here, g = J dtT is amplification gain, C 9 = cosh(g), T g = 
tanh(g). The notation |(2j) v , (2i + denotes 2j photons in polarization <p 

and 2i + 1 in ip- 1 , which contribute to the superposition with probability 
7 2 -. Fig. [2] depicts 7;o(*) 2 an d loj(j) 2 for g = 4 and reveals slow decay of 7 2 in 
infinity. Convergence of 7oj and 7^0 to zero is slower for higher gain. 

The states |$) and |<j>j_) are orthogonal due to different occupation parity 
in polarizations and are called macroscopic qubit. In experiment they con- 
tained 4m ~ 10 4 photons on average, where m = sinh 2 g. However, in this 
regime detection is not single photon resolving and they reveal effective overlap. 
Their distinguishability is efficiently quantified by photon number distributions 
p<s>(k,l) — |(fc,Z|<f>)| 2 and p$ x (k,l) giving probability of having k photons in 
polarization <p and I in (p 1 - simultaneously 

Mk, i) = p$_l (i, k) = H*-w/» for odd k and even l > ( 4) 

I otherwise, 

where X)fc°;=o P*(^> = 1- Distinguishability equals 

00 

v = l- ^2 vWM)p*l(M)- ( 5 ) 

k,l=0 

3.2. Additional operations performed on MQS 

Partial indistinguishability of |$) and is a drawback to any quantum 
protocol and Bell inequality test, and needs to be fixed by special quantum 
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state engineering. In 14, 1|| 17 1 a filtering method was suggested to increase 



distinguishability of macro-qubit. It is described by the projective measurement 

P< Tfc > = ]T \k,l){k,l\. (6) 

k,l;k+l>cr 

It cuts off the low photon number contributions below a threshold a in the 
initial superposition so that the preselected macro-states contain a photons at 
least distributed over two polarization modes. We will call this operation the 
theoretical preselection. Photon number distributions for MQS subjected to 
theoretical preselection read p^ h \k, I) = |(fc, l\T ,l £ : ' h \<&)\ 2 

n W(h l\ - Jl JVT ''l a Tfk-i)/2,i/2 for odd fc > even l,k + l>tr, , 
10 otherwise, 

P^\ k J) =P$> h \l-,k), normalization constant \N T h\~ 2 = YT+j>(<j-i)/2lij- 

For large class of useful observables, mean values evaluated for preselected 
multi-photon states will involve sums of the form | A^t^, | 2 X)°°_ i+j>a "fij /( 2 * + 
1, 2j), where f(p, q) is a polynomial. For example, for mean number of photons 
one sets f(p, q) = (p + q) and f(p, q) = (p + q) 2 for its variance. 

In the experiment, theoretical preselection is implemented by a weak mea- 
surement. The state is split into two beams, reflected and transmitted, by 
a beam splitter (BS) with high transmitivity. Quantum operation P^, is 
performed only on the reflected part and the result is feedforwarded to the 
transmitted beam (i.e. the summation constraint on occupations is shifted from 
reflected to transmitted part). This is the beam-splitter preselection. 

Action of a BS with reflectivity R (transmitivity T = 1 — R) on two orthog- 
onal polarizations ip and ip± is independent. For a given polarization, if one 
BS input is vacuum |0) and the other is photon number state \N), its action 
is given by Ubs\0,N) = J2k=o c i 1-^ — k) t \k) r with the probability amplitude 
cj^ = J {^)R k T N ~ k , where t (r) denotes mode where TV — K (k) photons 
were transmitted (reflected). For the multi-photon states we get Ubs\$) = 

E^=o7ij E?=o c " I+1) Em=o c ™ j) |( 2 * + 1 - n )v>( 2 j - ™)^)t\n vi m v ±.) r and 
similar expression for Ubs\®±)- The coefficients may be grouped as follows 

7ij c\ 2l+1 ^ = C 2 7io dh % ^ 7oj Cm^ ■ Evaluation of any physical quantity 
will lead to the following subexpressions 

MM) = c 2 ^ ( c f +d) 2 = cf (rp/af ^J^i-n)! ^ 1 ^ 1 "" < 8 > 

and fj{m,j) = C 2 j 2 j(cm ) 2 , where i,j,n,m are non-negative integers. For 
the beam-splitter preselection photon number distributions read 

oo 

P^ S \k,l) =\N B s\ 2 X] /j( m 'i) ' <5fc,2i+l-« Sl,2j-m, (9) 

i,j=0 0<n<2i+l 
0<m<2j 
n-\-m>(r' 
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with (k, I) = p^^Q, k), where normalization constant |A^ss| 2 = YlTj=o 

Y?n=o Sm=o,n+m>a' ' fj( m >j) and 6 a,b is a Kronecker delta function 

equal 1 if a = 6 and else. Here, important mean values are computed by 



\Nbs\ 2 EZ=o £n=o E^=o,„ +m >,' ' 3) f& + 1 - n, 2j - m). 

Beam splitter is useful for many other basic quantum operations, e.g. losses, 
homodyne detection, inefficient detection, entanglement distillation. 

3.3. Hypergeometric function 

We notice that the ratios (7i+i,o/7;,o) 2 an d fi(n,i + l)//,(n, i) are rational 
functions of i and i and n, respectively. According to the definition of a hy- 
pergeometric term [l^j ifj- fh fj arc double hypergeometric terms and their 
infinite sums over these parameters are hypergeometric functions. 

The sum of the probability governing MQS is given by Gaussian hypergeo- 
metric function 

G(n) =f> 2 = x n 2 F t Q^; z\ , G{m) 7* = Dm 2 f i (]^ v z) , (10) 

where x n = (2n + l)\T 2n /{C 2 n\ 2 2 2 ™), y m = (2m)! T 2m /(C 2 g ml 2 2 2m ), z = T 2 , 
a = n + | , 6 = m+ i. The situation is similar for the probability governing the 
MQS after passing BS. Here, the sums A(n) = Y^iL\n/ 2 \fi(. n ^) and B(m) = 
Y^j=\(m+l)/2\fj( m ->j)-> wncrc l x \ is the floor function [a; J = max{m € Z|m < 
x}, have to be written for odd and even n and m separately 



A(2n) =x n (2n+l)R 2n T 2 F 1 ( a f;T 2 z) , A(2n + 1) = x n R 2n+1 2 F X T;T 2 



" z 



2 

(11) 



B(2m) =y m R 2m 2 F 1 ^T;T'zj ,B(2m- 1) = y^mR^T 2 F t (V;!"* 

We checked that no closed form of G(n), G(m), A(n) and B(m) exist using 
the Gosper's, Zeilberger's, and Petkovsek's algorithms [l8|. These algorithms 
constitute the standard mathematical tools for hypergeometric series analysis. 
Closed form, if existed, would eliminate sums over i and j completely and reduce 
the computation time. 

In our case, these sums have to be computed itcratively. This prevents effi- 
cient simulation of MQS properties and evolution for realistic parameters, since 
it is neither possible to bring this function to elementary ones nor to other 
special functions (e.g. Bessels) tractable more easily. The function 2 F\ arises 
directly form probability governing the MQS 7^. Parameters a and b of 2 Fi are 
half integers. This class of Gaussian hypergeometric functions is the least known 
in the literature. There are only a few known identities simplifying 2 Fi with 
half- integer parameters [3, Gil, but they cannot be used here. Moreover, re- 
currence transformations lead to other half-integer forms of 2 Fi. Unconditional 
convergence of the hypergeometric sum is guaranteed by z, T 2 £ (0, 1). 
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It is difficult to find the appropriate cutoffs for the sums over i and j. High 
cutoffs make the computation long and error prone whereas too low values result 
in omitting significant terms. The fast computation of 2-F1 in a non-asymptotic 
regime of its arguments and parameters is challenging 2(| 21 1. In literature, 



there are al gori thms applicable for some special cases of relations between its 
parameters (l{§ |22|, |23|, but not useful here. Moreover, in our intermediate 



regime they become unstable. Neither Gaussian quadrature [24| nor Pade [25 
approximations exist. 

To be able to compute Eq. © and similar ones of the form 

00 2i+l 2j 

5=EE E fMfMJ) (12) 

i,j= n— m—0;n-\-m>cr' 

we aim at partial factorization with respect to (i,n) and (j,m) to compute 
the sums over i and j separately. Full factorization cannot be achieved due 
to the preselection condition n + m > a'. We change the variable order 
E^oErf/.fn^) = Er=oE£[n/2J M n ^) (similarly for fj(m,j)) inEq. CE 



A{n)B{m). (13) 



n,m— 0;n+m>cr' 

We avoid four nested sums and get a double sum of independent series. How- 
ever, it became clear, that the hypcrgcomctric functions have to be additionally 
summed up in infinite sums over its parameters. 

In case of computation of mean values, f(p, q) has first to be split into 
sum of monomials f(j>, q) = J2 n m a nmP n 9™, where, importantly, a nm are real 
coefficients. Then hypergeometric terms take the form /i(n, i) = fi(n, i) ■ (2i + 
1 — n) p and fj(m 7 j) = fj(m,j) ■ (2j — m) q , where p and q are integers. The 
analysis presented above for fi{n,j) still holds for /,; and fj. 



4. Numerical Methods 

4-1- Calculation of Hypergeometric Terms with High Precision 

Factorials in numerators and denumerators of hypergeometric terms 7^0 , joj , 
fi and fj take large values, unavailable for the standard machine arithmetic. 
IEEE 754 double precision numbers are capable of storing 15 significant decimal 
digits and magnitudes 1CT 308 to 10 308 0. For i, j > 50 they exceed 10 308 , but 
the values of terms are between and 1 and need to be computed precisely. 

One way to compute jf , 7^, fi and fj is to use libraries offering enhanced 
precision arithmetic. For example, the GNU Multiple Precision Arithmetic 
Library (GMP) l27ll and the GNU Multiple Precision Floating-Point Reliable 
Library (MPFR) 12811 are available for the C and C++. The Class Library for 
Numbers (CLN) }29| ] cooperates with C++ and the mpmath [3(| with Python. 
These libraries offer low speed compared to the double precision arithmetic. 
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Figure 3: Function fi(n,i) evaluated for i, parametrized by n (g = 4). 



The other possible way of computing hypergeometric terms is to use their 
logarithmic forms: log %q, log7(y, log/i, log/j. The number magnitudes are 
within range of double precision for large i and j. Moreover, both the standard 
C/C++ and MPFR offer lgamma(n) function, which calculates log(r(n)) = 
log((n — 1)!) with hardware acceleration. 

In order to find the best method of computation of the subexpressions we 
compared all the mentioned solutions. The relative error estimation was based 
on the values obtained symbolically with computer algebra system Wolfram 
Mathcmatica. We showed superiority of C/C++ computation with double- 
precision. It is not only several orders of magnitude faster but also consumes less 
memory. The reason for that is that IEEE 754 arithmetic is built-in into modern 
computer processors and does not need any additional memory structures. The 
average error of double-precision calculation is of the order of 10 -12 , low enough 
for numerical simulations. The error can be further decreased at the significant 
cost of speed with MPFR library, which is the second fastest solution. 



4-2. Convergence Rate 

Computation of G(n), G(m), A{n) and B{m) requires finding proper cutoffs 
for infinite sums over i and j in Eqs. (|10[) and First, we tested con- 

vergences of these sums, since their acceleration allows setting cutoffs earlier. 
Several convergence rate acceleration methods were tested: the Aitken's delta- 
squared process 3ll 321. the Shanks transformation [33| and the Richardson 
extrapolation 34j, [35j. None of them improves much the computation time. 



They raise the complexity of the algorithm instead. 



4-3. Efficient Computation of the Gaussian Hypergeometric Function 

Next, we worked out algorithm for finding cutoffs for sums A(n) and B{m), 
but simplified procedure holds for G(n) and G(m). The hypergeometric terms 
7? , 7q-, fi(n,i) and fj(m,j) converge slowly. Function fa is depicted in Fig. [3] 
(see Fig. [2] for jij). The greater n gets, the slower /j converge over variable 
i and more terms have to be included in the summation (similarly for fj and 
variable j). One may apply high constant sum cutoffs. But even for large i (j) 
there exist n (m) for which the fixed range is too small. 

As a solution to this problem, we applied a dynamic cutoff with an adjustable 
precision e.g. 10 -15 . The summation range was divided into two intervals. In the 
first one, the terms fi and fj monotonically grow to achieve the global maximum. 
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Table 1: Comparison of cutoff values and relative error of computation of A(n) for g = 4 with 
the fixed and dynamic cutoff estimation methods. 

The second interval is infinite and the terms slowly decrease asymptotically to 
0. The value i = I for which /j achieves its maximum can be found by solving 
differential equation a ^g"'^ = 0, which leads to equation ^q(2I + 1 — n + 1) = 
^o( n + 1) + log R ~ l°g T, where ^o( x ) denotes the digamma function (the first 
derivative of logr(x)). The inverse of ^0(2/ + 1 — n + 1) can be computed 
numerically with the Paul Fackler's method (36j . Similar method may be used 
for finding maximum of fj. 

Alternatively, algorithm starts from i = \ n/2\ (j = L(m+1)/2J) and for each 
i (j) evaluates fi{n,i) (fj(m,j)). It adds this term to the sum and compares 
fi(n,i) to the previous term fi(n,i — 1). If fi(n,i) > fi(n,i — 1) the algorithm 
is still in the first interval and continues summation. Otherwise, the summation 
is performed over the second interval, the algorithm checks if fi is smaller than 
the desired precision and if yes, it stops. Error is minimized by adding the terms 
in the increasing order in the first interval. In practice, the algorithm uses the 
logarithmic forms (they allow finer comparison of fi(n,i) and fi(n,i — 1)). 

We compared the fixed and dynamic methods of assigning cutoffs for the 
sum A(n) in Table [T] (the numbers for B(m) are similar). The latter decreases 
the number of terms for small n and m and achieves a good precision for their 
large values. It significantly improves speed and accuracy of the computation 
e.g. for g = 4, if n = 100, the summation includes 1000 whereas if n = it 
includes 150 terms to get the same accuracy with relative error 10~ 15 . 

fi and fj have the same properties as fi and fj, only converge slower. There- 
fore, the presented algorithm is capable of finding cutoffs of infinite hypergeo- 
metric sums required in computation of mean values. 
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4-4- Summation of Hypergeometric Functions 

Computation of the sum S in Eq. (|13p involves infinite summation of product 
of two Gaussian hypergeometric functions A(n) and B(m) over n and m, where 
n + m > a' . The solution is to pre-compute the values of A(n) and B(m) with 
method presented in Section 14.31 and store them in the computer memory. The 
next step is to perform the final summation. The cutoffs TV and M over variables 
n and m, respectively, are found separately during the pre-computation stage, 
so that \A(N)/A(a')\ and \B(M) / B(a')\ are smaller than the desired relative 
computation error. 

Additional speedup of computations is achieved by changing the order of 
summation and eliminating one of the infinite sums 

oo n 

S=Y^ ^2A(n~m)-B(m). (14) 

n—a' m— 

The computation time of this formula is polynomial 0(N 2 ). 

Further optimization is possible by noting the fact that hypergeometric 
terms jf j: f 4 , fj are probability distributions and Yli°j=o 7ij = ^ J2ij=o J2l=o 

Em=o fi( n ' /i( m '.?) = !> wnicn implies X)~ m=0 A ( n ) B ( m ) = 1 - Thcsc P r0 P" 
erties allow to rewrite Eq. (|14[) into a form involving finite summation range 

<T a —n 

5=l-^A(n)^B(m). (15) 

n— m— 

Here computation time is also polynomial 0(o~' 2 ). 

Eq. (|15p takes advantage over Eq. (fT4]) for small values of a', for which the 
sum intervals are relatively short and the result is obtained quickly. However, 
the larger value of a' is, the more terms in Eq. (fT5|) have to be taken into 
account. As a result of that, the summation takes longer and errors accumulate 
significantly. At some point Eq. (|14p gives more reliable results and performs 
faster than Eq. (|15|) . Thus, the routine computing S decides, depending on a', 
which method to apply. 

5. Examples and Applications 

5.1. Computation of MQS Normalization and Photon Number Distribution 

Formulae for photon number distributions, e.g. Eq. ©or ([§]) require prior 
computation of normalization |Aj^| 2 and \Nbs\ 2 , respectively. In principle, 
it is possible to take advantage of the fact that J2'k'i=oP^(^' = compute 
unnormalized distribution values for all pairs (k,l), sum them up to obtain 
the inverse of its normalization constant and later renormalize the result. This 
method is useful in cases where all meaningful values of p$(k,l) have to be 
computed anyway. However, in general, the normalization constant should be 
computed separately. 
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For theoretical preselection, according to Eqs. (JT3J) and (fT5)) , normalization 
constants equal 

oo i 

C tH £^-i),o"7o 2 ;, (16) 

i=(<r-l)/2 i=o 

(cr-l)/2 (<r-l)/2-i 

E E 7oV (17) 

t=0 3=0 

We empirically found for g = 4, a w 5000 to be a good threshold for switching 
from Eq. (|T7|) to Eq. (|T6|) . According to Eq. (J7J, computation of |-/Vr/i| 2 directly 
gives p$\k,l). 

\Nbs\ 2 is computed directly with Eqs. (fl"4"]l and ([15]) , We found that for 
.9 = 4 and a' w 500 the computation time required by Eqs. (|T4"|) and (fl"5)) is 
similar and for tr' > 500 the former is faster. In order to compute photon 
number distribution p^ S \k, I) we turn Eq. (|9]) into a simpler form 

oo n [ fi ( fc+n ~ m ~ 1 ,n-m) if k + n - m is odd 
pi BS) (k,l) = \N BS \ 2 EES -fi (^'™) and Z + m is even, 

n=CT ' m=o otherwise. 

(18) 

It is similar to Eq. (fl~6)) and is evaluated likewise. The main difference is that 
sums in Eq. (|18p traverse only even or odd values of n and m, depending on the 
parity of k and I. The computation time is polynomial 0(N 2 ), where N is the 
cutoff of the infinite sum over n. 



\N T h\ A 



\N Th \ 



5.2. Computation of MQS Distinguishability 

Distinguishability given in Eq. ([S]) is not optimal for computation, since it 
does not minimize summation errors and for large values of the cutoffs it is 
difficult to keep coefficients of the photon number distributions in memory. 

5.2.1. Optimization of Mathematical Formulae 
Eq. ([5]) may be optimized by change of variables 

K fc-l K 

v = i - 2 E VMfc> o p*q, *) - E^( fc > ( 19 ) 

k=0 1=0 k=0 

where K is the cutoff of the infinite sum over variable k. The computation time 
of Eq. ([T^| is polynomial 0(K 2 ) and has the same advantages and disadvantages 
as Eq. ([5]), but requires only half of the calculations. 

More complex idea is to first, divide the (k, I) plane into rectangles, compute 
the photon number distributions for each rectangle separately and save the data 
to a disk. Next, we calculate partial overlaps for each rectangle and sum up the 
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Figure 4: The plot of p^ SS '(fc,0) showing three main intervals of k where photon number 
distribution behaves differently. 

contributions [37]] • Our approach uses tiles of a size K' x V and performs 
two-step summation 

K/K' L/L' /K'-l L'-l \ 

v=i- E EE Ev^m»mL ( 20 ) 

fc=0 1=0 \k'=0 l'=0 ) 

where x = K'k + k', y = L'l + V and K , L are the cutoffs of the sums over k and 
I respectively. Computation complexity of Eq. (f2"0")) is the same as of Eq. ([5]), but 
this method offers several advantages. The values of K 1 and L' can be selected in 
such a way that the precomputcd values of coefficients (e.g. A(n) and B{m)) are 
stored in memory and used for all points of the single tile. Additionally, different 
tiles can be computed in parallel with separate processors. This significantly 
speeds up the calculations and leads to a better usage of the computer resources 
and minimizes errors by grouping summed terms in smaller sets, but requires 
proper estimation of cutoffs K, L and tile size K' x L'. Computation gets twice 
shorter if Eq. ([2"U]) is rewritten to the following form 



K/K' fc-l /K'-l K'-l \ 

11=1-2 E E E E \Jp<s>{x, w )p<s>{w,x) ) 

fc=o ;=o \ k'=o i'=o / 

K/K' /K'-l k'-l 

- E ( E x ) + 2 E \fpi>{ x , w )pt>{ w , x ) 

k=0 \k'=0 l'=0 



(21) 



where w = K'l + I'. 



5.2.2. Approximation: Cutoffs of Infinite Sums 

In order to estimate cutoffs K and L in Eqs. (|20| and (f2~Tj) we tried the 
method presented in Section T4.3I However, unlike A(n) and B(m), Eq. (|20[) 
involves a lot of near-zero terms for small k and I before entering the more 
significant regions and we had to adjust this algorithm accordingly. The new 
solution analyzes the shape of photon number distributions. The distribution 
cut for a given I spans of three different intervals of k, see Fig.UJ The first inter- 
val (A) starts at k = and ends at k rj 10 a' for the beam-splitter preselection 
(10 a for the theoretical). Here, the values of the distributions are small and 
their input is negligible. The interval B for k £ (10er',30m) (k £ (10 cr, 20m)) 
gives the main input. In region C the values of are slowly fading out in infinity. 

The proposed algorithm analyzes the values of p$(fc,0) to find the start of 
the region C. To speed up the search, it begins at k = 30m (k = 20m), which 
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is an approximate end of region B, and takes samples of the photon number 
distribution for k = 30m, 30m + 100, 30m + 200 ... to the point where the values 
of the distribution are not distinguishable from with the desired precision. 
This point lays inside region C and can be used as a good cutoff value K and L. 
This method gives the results which confirm the values obtained analytically. 

5. 3. Simulation of Finite Detector Resolution 

For the theoretical preselection we model the resolution of detectors mea- 
suring photon number ±150 with the Weierstrass transform, i.e. we apply the 
low-pass filter with Gaussian distribution for which 3cj ss 150. The photon 
number distribution p'^ h \k, I) takes the form 



where a is the standard deviation. Gaussian properties of the filter allow for 
replacing the infinite summation over p and q with finite ranges. 

In case of the beam-splitter preselection, applying a Gaussian filter to the 
photon number distribution values does not influence the result in a signif- 
icant way. Since the filtering slows the computations, which in case of the 
beam-splitter preselection are already time-consuming, there is no reason for 
the Weierstrass transform. 

Application of Eq. (|2"2"|) to Eq. (|2"Tj) means that each tile of size K' x K' has 
to be increased by a 3ct margin. Computation time for each tile is increased 
approximately (1 + 6a /K') 2 times and requires more computer memory, but is 
still quite easily tractable. 

6. MQSVIS Program Suite 

MQSVIS is a program suite developed for computation of various quantum 
indicators for MQS of light. These indicators include: 

• photon number distributions for MQS of light, 

• visibility, 

• visibility computed from the overlap, 

• visibility computed for the photon number distributions transformed with 
Weierstrass transform for different values of 3<r 6 {1, 1.5, 15, 150} (Gaus- 
sian blur), 

• mean number of photons in both polarizations (jointly and separately for 
each polarization), 

• variance of number of photons in both polarizations (jointly and separately 
for each polarization). 




(22) 
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6.1. Compilation of the Programs 

In order to compile and test MQSVIS program suite for Linux/Unix op- 
erating systems with GNU utilities (GNU make, GNU compiler collection) it 
is enough to run the following commands in a directory containing unpacked 
source code 

make 

make check 

In case of Linux/Unix with GNU make and Intel C Compiler, one has to modify 
Makefile, replacing gcc with ice and -f openmp option with -openmp. 

For other platforms and compilers: one should compile mqsvis_norm. c with 
a standard C compiler, and mqsvis_tile . c with a C compiler and OpenMP ex- 
tensions turned on. If compiler docs not offer OpenMP extensions, the program 
will still work but will not utilize multiple cores or processors. 

6.2. Description of the Program Suite 

6.2.1. MQSVIS_norm 

Program computes squared normalization for preselected macroscopic quan- 
tum superpositions (MQS) of light. This value is required by MQSVISMle 
program. 

Program parameters: m - average number of photons, Dth - preselection 
threshold. 

6.2.2. MQSVISMle 

Program computes photon number distribution for preselected macroscopic 
quantum superpositions (MQS) of light for a single tile. It produces partial 
results of quantum indicators: visibility, man number of photons, variance, 
to be gathered by MQSVIS-gather script. Additionally, saves photon number 
distribution for single MQS state to a hie. 

Program parameters: m - average number of photons, Dth - preselection 
threshold, R - amount of losses, N2 - squared normalization, computed by 
MQSVIS-norm program, tilesize - width and height of a single tile, tilex - tile 
column (counted from 0), tiley - tile row (counted from 0), plotstep - step used 
for saving photon number distributions plotfnamel - file to save photon num- 
ber distribution for tile (tilex, tiley), plotfname2 - file to save photon number 
distribution for tile (tiley, tilex). 

Environment variables (used only when compiled with OpenMP extensions): 
OMP_NUM_THREADS - maximal number of cores or processors utilized by the 
program. 

6.2.3. MQSVIS.gather.py (Python script) 

A script developed to gather partial results computed by MQSVISMle pro- 
gram. It computes final quantum indicators for macroscopic quantum superpo- 
sitions of light. 

Program requires partial results of MQSVISMle program to be saved in text 
files with the names in the form 
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M°/„s_Dth /„s_r'/.s- /.d , 7.d . txt 



where first %s is replaced with average number of photons, second °/„s with 
preselection threshold, third °/„s with losses and °/ d, '/A arc comma-separated tile 
coordinates numbered from 0. 

Program parameters: m - average number of photons, Dth - preselection 
threshold, R — losses, tiles - number of tiles in a row/column. 

6.3. Examples and Testing 

Below we present a sample session with the above programs. The sequence 
of commands could be used e.g. for testing of the suite. 

1. Compute normalization for m = 5 and Dth = 2 and save it to a file 

mqsvis_norm 5 2 > normalization.txt 

2. Compute tile (0,0) for R = (no losses), tile size 10 x 10 

mqsvis_tile 5 2 $(< normalization.txt) 10 1 \ 
plot-O.O.txt /dev/null > M5_Dth2_rO-0,0.txt 

3. Compute tile (1,0), and parallely tile (0, 1) for the same parameters 

mqsvis_tile 5 2 $(< normalization.txt) 10 1 1 \ 
plot-l.O.txt plot-0,l.txt > M5_Dth2_rO-l,0.txt 

4. Compute tile (1, 1) 

mqsvis_tile 5 2 $(< normalization.txt) 10 1 1 1 \ 
plot-l.l.txt /dev/null > M5_Dth2_rO-l , 1 .txt 

5. Gather the results 

python mqsvis_gather .py 5 2 2 

The computed results for the above session, printed out in the last step, are 
as follows 

total probability sum=0 . 626948016783147 
visibility=0 . 683 
visibility prim=0.683 

simple visibility computed from the overlap=l . 000 

visibility with Gaussian blur (3sigma=l) =0.968 

visibility with Gaussian blur (3sigma=l . 5)=0 . 746 

visibility with Gaussian blur (3sigma=15)=0.491 

visibility with Gaussian blur (3sigma=150)=0.912 

mean=8 . 362 

mean k=6.269 

mean 1=2.092 

variance=68 . 303 

variance k=40.241 

variance 1=15.715 

maximal value=0 . 0408525598455265 
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7. Conclusions 



We discussed numerical methods useful for modeling of macroscopic quan- 
tum superpositions of light generated by phase covariant quantum doner. Prop- 
erties of these MQS are governed by Gaussian hypergeometric function, which 
parameters are half integers and argument takes moderate values. No simplifica- 
tions of this function are possible. It is given by slowly convergent infinite sum, 
impossible to accelerate. The algorithm dynamically finds cutoff as well as cut- 
offs of sums involving hypergeometric functions. We achieved precision 10~ 15 . 

The model allows simulation of experimental schemes involving linear opti- 
cal elements and MQS, as well as computation of expectation values of some 
observablcs. It is not possible to obtain such results by analytical calculation. 
Our model takes into account values of experimental parameters. Depending 
on them optimizes computation by choice of the more efficient, parallelizablc 
algorithm for evaluation of the indicators for these states, e.g. the photon num- 
ber distribution or distinguishability. Numerical methods presented here can be 
adjusted for analysis of decoherence and similar experimental schemes, different 
filtering techniques and modified MQSs. As an example, we included a model 
of a realistic detector what improved applicability of our model for scientific 
research. We verified our numerical results with recently developed analytical 
model for MQS, available only for certain regimes of parameter values [14j. 




Finally, MQSVIS program suite was developed and implemented in C and 
Python programming languages. It utilizes the developed numerical model and 
allows to compute various quantum indicators for MQS of light: photon number 
distributions, visibility in case of limited resolution of photodctcctors, mean 
and variance of photon number (jointly and separately for each polarization). 
Parallelization was achieved with splitting the computing tasks into tiles and 
utilization of OpenMP compiler extensions. 
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